Dataset Exploratory Analysis

Module-Trait associations

datTraits = datMeta %>% dplyr::select(Diagnosis_, Region, Sex, Age, PMI, RNAExtractionBatch) %>%
            rename('Diagnosis' = Diagnosis_, 'ExtractionBatch' = RNAExtractionBatch)

# Recalculate MEs with color labels
ME_object = datExpr %>% t %>% moduleEigengenes(colors = dataset$Module)
MEs = orderMEs(ME_object$eigengenes)

# Calculate correlation between eigengenes and the traits and their p-values
moduleTraitCor = MEs %>% apply(2, function(x) hetcor(x, datTraits, ML=TRUE)$correlations[1,-1]) %>% t
rownames(moduleTraitCor) = colnames(MEs)
colnames(moduleTraitCor) = colnames(datTraits)
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nrow(datExpr))

# In case there are any NAs
if(sum(!complete.cases(moduleTraitCor))>0){
  print(paste0(sum(is.na(moduleTraitCor)),' correlation(s) could not be calculated')) 
}

## [1] "1 correlation(s) could not be calculated"

moduleTraitCor = moduleTraitCor[complete.cases(moduleTraitCor),]
moduleTraitPvalue = moduleTraitPvalue[complete.cases(moduleTraitCor),]

# Sort moduleTraitCor by Diagnosis
moduleTraitCor = moduleTraitCor[order(moduleTraitCor[,1], decreasing=TRUE),]
moduleTraitPvalue = moduleTraitPvalue[order(moduleTraitCor[,1], decreasing=TRUE),]

# Create text matrix for the Heatmap
textMatrix = paste0(signif(moduleTraitCor, 2), ' (', signif(moduleTraitPvalue, 1), ')')
dim(textMatrix) = dim(moduleTraitCor)


labeledHeatmap(Matrix = moduleTraitCor, xLabels = names(datTraits), yLabels =  gsub('ME','',rownames(moduleTraitCor)), 
               yColorWidth=0, colors = brewer.pal(11,'PiYG'), bg.lab.y = gsub('ME','',rownames(moduleTraitCor)),
               textMatrix = textMatrix, setStdMargins = FALSE, cex.text = 0.8, cex.lab.y = 0.75, zlim = c(-1,1),
               main = paste('Module-Trait relationships'))

rm(ME_object, datTraits, MEs, moduleTraitCor, moduleTraitPvalue, textMatrix)

Selecting the three modules with highest (absolute) correlation to Diagnosis: #D675FD, #00BF78 and #13B700

pca = datExpr %>% prcomp

plot_data = data.frame('ID'=rownames(datExpr), 'PC1' = pca$x[,1], 'PC2' = pca$x[,2]) %>%
            left_join(dataset, by='ID') %>% dplyr::select(ID, PC1, PC2, Module, gene.score) %>%
            mutate(ImportantModules = ifelse(Module=='#D675FD', '#D675FD',
                                      ifelse(Module=='#00BF78', '#00BF78',
                                      ifelse(Module=='#13B700', '#13B700', 'Others')))) %>%
            mutate(color = ifelse(ImportantModules=='Others','gray',ImportantModules),
                   alpha = ifelse(ImportantModules=='Others', 0.1, 0.5))

table(plot_data$ImportantModules)

## 
## #00BF78 #13B700 #D675FD  Others 
##    1158     182    1277   13983

ggplotly(plot_data %>% ggplot(aes(PC1, PC2, color=ImportantModules)) + 
         geom_point(alpha=plot_data$alpha, color=plot_data$color) + theme_minimal() + 
           ggtitle('Modules with strongest relation to Diagnosis'))

Gene Significance

Comparing this with the DEA plot, Gene significance close to 0 means the gene expression is not affected by Diagnosis and negative values means the gene is underexpressed in Autism samples. So the important thing about gene significance is its absolute value

plot_data = data.frame('ID'=rownames(datExpr), 'PC1' = pca$x[,1], 'PC2' = pca$x[,2]) %>%
            left_join(dataset, by='ID')

ggplotly(plot_data %>% ggplot(aes(PC1, PC2, color=GS)) + geom_point(alpha=0.5) + 
         scale_color_gradientn(colours=c('#F8766D','white','#00BFC4')) +
  theme_minimal() + ggtitle('Gene Significance'))

Module Membership and Gene Significance relation

plot_data = dataset %>% dplyr::select(MM.D675FD, GS, gene.score) %>% filter(dataset$Module=='#D675FD')

ggplotly(plot_data %>% ggplot(aes(MM.D675FD, GS, color=gene.score)) + geom_point(alpha=0.5) + 
         scale_color_manual(values=SFARI_colour_hue(r=c(1:6,8,7))) + 
         theme_minimal())

plot_data = dataset %>% dplyr::select(MM.00BF78, GS, gene.score) %>% filter(dataset$Module=='#00BF78')

ggplotly(plot_data %>% ggplot(aes(MM.00BF78, GS, color=gene.score)) + geom_point(alpha=0.5) + 
         scale_color_manual(values=SFARI_colour_hue(r=c(1:6,8,7))) + 
         theme_minimal())

plot_data = dataset %>% dplyr::select(MM.13B700, GS, gene.score) %>% filter(dataset$Module=='#13B700')


ggplotly(plot_data %>% ggplot(aes(MM.13B700, GS, color=gene.score)) + geom_point(alpha=0.5) + 
         scale_color_manual(values=SFARI_colour_hue(r=c(2:6,8,7))) + 
         theme_minimal())

Identifying important genes

Selecting the modules with the highest correlation to Diagnosis, and, from them, the genes with the highest module membership-(absolute) gene significance

*Ordered by \(\frac{MM+GS}{2}\)

Only three SFARI Genes in the three lists …

top_genes_1 = dataset %>% dplyr::select(ID, MM.D675FD, GS, gene.score) %>% filter(dataset$Module=='#D675FD') %>%
              rename('MM' = 'MM.D675FD') %>% mutate(importance = (MM+abs(GS))/2) %>% arrange(by=-importance) %>%
              top_n(10)

kable(top_genes_1, caption='Top 10 genes for module #D675FD')

top_genes_2 = dataset %>% dplyr::select(ID, MM.00BF78, GS, gene.score) %>% filter(dataset$Module=='#00BF78') %>%
              rename('MM' = 'MM.00BF78') %>% mutate(importance = (MM+abs(GS))/2) %>% arrange(by=-importance) %>%
              top_n(10)

kable(top_genes_2, caption='Top 10 genes for module #00BF78')

top_genes_3 = dataset %>% dplyr::select(ID, MM.13B700, GS, gene.score) %>% filter(dataset$Module=='#13B700') %>%
              rename('MM' = 'MM.13B700') %>% mutate(importance = (MM+abs(GS))/2) %>% arrange(by=-importance) %>%
              top_n(10)

kable(top_genes_3, caption='Top 10 genes for module #13B700')

plot_data = data.frame('ID'=rownames(datExpr), 'PC1' = pca$x[,1], 'PC2' = pca$x[,2]) %>%
            left_join(dataset, by='ID') %>% dplyr::select(ID, PC1, PC2, Module, gene.score) %>%
            mutate(color = ifelse(Module=='#D675FD', '#D675FD',
                           ifelse(Module=='#00BF78', '#00BF78',
                           ifelse(Module=='#13B700', '#13B700', 'gray')))) %>%
            mutate(alpha = ifelse(color %in% c('#D675FD', '#00BF78', '#13B700') & 
                                  ID %in% c(as.character(top_genes_1$ID), 
                                            as.character(top_genes_2$ID),
                                            as.character(top_genes_3$ID)), 1, 0.1))

plot_data %>% ggplot(aes(PC1, PC2)) + geom_point(alpha=plot_data$alpha, color=plot_data$color) + 
              theme_minimal() + ggtitle('Important genes identified through WGCNA')

SFARI genes and Module-Diagnosis correlation

Modules with the strongest module-diagnosis correlation should have the highest percentage of SFARI Genes, but the plot shows this relation in modules with correlations close to 0, and the behaviour is the same for the other clusterings I have tried…

A good thing is that the gray cluster (which is actually all the genes that were left without a cluster) has a low correlation with ASD and also a low percentage of SFARI genes

plot_data = dataset %>% mutate('hasSFARIscore' = gene.score!='None') %>% 
            group_by(Module, MTcor, hasSFARIscore) %>% summarise(p=n()) %>% 
            left_join(dataset %>% group_by(Module) %>% summarise(n=n()), by='Module') %>% 
            mutate(p=round(p/n*100,2)) 

for(i in 1:nrow(plot_data)){
  this_row = plot_data[i,]
  if(this_row$hasSFARIscore==FALSE & this_row$p==100){
    new_row = this_row
    new_row$hasSFARIscore = TRUE
    new_row$p = 0
    plot_data = plot_data %>% rbind(new_row)
  }
}

plot_data = plot_data %>% filter(hasSFARIscore==TRUE)

ggplotly(plot_data %>% ggplot(aes(MTcor, p, size=n)) + geom_smooth(color='gray', se=FALSE) +
         geom_point(color=plot_data$Module, aes(id=Module)) +
         xlab('Module-Diagnosis correlation') + ylab('% of SFARI genes') +
         theme_minimal() + theme(legend.position = 'none'))

rm(i, this_row, new_row)